Quantum Monte Carlo and variational approaches to the Holstein model 
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Based on the canonical Lang-Firsov transformation of the Hamiltonian we develop a very effi- 
cient quantum Monte Carlo algorithm for the Holstein model with one electron. Separation of the 
fermionic degrees of freedom by a reweighting of the probability distribution leads to a dramatic 
reduction in computational effort. A principal component representation of the phonon degrees 
of freedom allows to sample completely uncorrelated phonon configurations. The combination of 
these elements enables us to perform efficient simulations for a wide range of temperature, phonon 
frequency and electron-phonon coupling on clusters large enough to avoid finite-size effects. The 
algorithm is tested in one dimension and the data are compared with exact-diagonalization results 
and with existing work. Moreover, the ideas presented here can also be applied to the many-electron 
case. In the one-electron case considered here, the physics of the Holstein model can be described 
by a simple variational approach. 

PACS numbers: 63.20.Kr, 71.27.+a, 71.38.-k, 02.70.Ss 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) simulations for mod- 
els with electron-phonon coupling are often limited in 
both system size and accessible parameter range by long 
autocorrelation times and large statistical errors. This 
makes it very difficult to study realistic models for, e.g., 
the high temperature superconductors or the manganites 
which exhibit colossal magnetoresistance. In both classes 
of materials electron-phonon interactions play an impor- 
tant rolei^ Although classical treatments of phonons 
have been quite successful in certain situations, 3,4 e.g., at 
high-enough temperatures, quantum effects are expected 
to be relevant^ Consequently, it is highly desirable to de- 
velop a new, more efficient method to treat the phonon 
degrees of freedom quantum mechanically. The Holstein 
molecular-crystal model constitutes one of the simplest 
models for coupled electron-phonon systems, and there- 
fore serves as an ideal testing ground for new approaches. 
Moreover, despite enormous theoretical efforts, even the 
physics of the Holstein model is still not completely un- 
derstood. 

The extensive use of QMC methods to study many- 
body problems is based on the fact that they can give 
quasiexact results (i.e., exact apart from statistical er- 
rors which can be made arbitrarily small, in principle). 
Over the years, several different QMC methods have been 
developed to study systems with electron-phonon cou- 
pling, such as the Holstein^ the Fr6hlich)& or the Su- 
Schrieffer-Heeger (SSH) modeli A very general QMC 
method for coupled fermion-boson models was developed 
by Blankenbecler et al& and Scalapino and Sugar It is 
based on an analytic integration over the fermion de- 
grees of freedom and a MC simulation of the resulting 
boson model. The simulation is performed using the 
grand-canonical ensemble and requires the evaluation of a 
fermion determinant involving a computation time which 
is proportional to the cube of the system size. Moreover, 
the method in its original form becomes unstable at low 
temperatures. While the simulations of Refs. H and 



were restricted to one dimension, Levine and Su-**** and, 
using a stabilized version of the same algorithm, Niyaz 
et alM> studied charge-density-wave formation and su- 
perconductivity in the two-dimensional Holstein model. 
A numerically faster method is the world-line algorithm 
developed by Hirsch et aZ^2»i» based on a special breakup 
of the Hamiltonian and a fixed number of fermions. This 
results in configuration weights which are simple to eval- 
uate allowing for much bigger system sizes. In the course 
of the simulation, both fermions and bosons are sam- 
pled simultaneously. The latter method has been suc- 
cessfully applied to the Holstein polaron problem and to 
the half-filled SSH and Holstein modeli 13 i 14 i 15 i 16 i 17 How- 
ever the world-line algorithm is restricted to models in 
one spatial dimension or to the single-electron case in 
any dimension by the minus-sign problem. 18 Scalettar 
et alm^ applied a rather complicated so-called hybrid 
molecular dynamics algorithm to the two-dimensional 
Holstein model near half filling. This work was ex- 
tended to the low-temperature regime by Noack et alm^ 
Finally, MarsiglioSi developed a low-temperature QMC 
method to study the same model, also at half filling. De 
Raedt and Lagendijk22i2i24 and Kornilovitch 2 ^ used an 
alternative approach based on Feynman's path-integral 
method^ where the boson degrees of freedom are inte- 
grated out analytically and the resulting fermionic model 
is simulated using QMC. Although the method is limited 
to one electron or two electrons of opposite spinal by the 
sign problem, it allows efficient simulations in one, two, 
and three dimensions even for small phonon frequencies 
near the adiabatic limit, and has also been used to inves- 
tigate the Holstein model with dispersive phonons. 24 Also 
using Feynman's path integral, Kornilovitch and Pike 28 
developed a QMC method which exploits the conserva- 
tion of the total quasimomentum of the system and al- 
lows the calculation of dynamical properties such as, e.g., 
the polaron band structure. Although the method is not 
restricted to a certain model or dimensionality of the sys- 
tem, it suffers from large statistical errors. Moreover, it 
is limited to the case of a single fermion at very low tern- 
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perature, and exhibits a sign problem for nonzero total 
quasimomentum even in one dimension. Prokof'ev and 
Svistuno\«2& and Mischenko et al£^ used QMC to directly 
sample the zero-temperature one-electron Green function 
of the Frohlich polaron in imaginary time. The method 
allows calculations for an infinite system in three dimen- 
sions, but requires a convergent series for the electron 
propagator. While all but the last method mentioned 
so far make use of the Trotter-Suzuki approximation, 18 
Kornilovitcb^i^ developed a continuous-time algorithm 
that works in any dimension and allows calculations on 
infinite systems. It gives directly dynamical quantities 
such as the polaron spectrum and effective mass with 
very high accuracy. However, similar to previous work^S 
it is restricted to the one-electron limit at very low tem- 
peratures. Moreover, calculations for small phonon fre- 
quencies and/or weak electron-phonon coupling are diffi- 
cult and a sign problem appears for nonzero total quasi- 
momentum. The projector QMC method 18 in combi- 
nation with a local updating of the phonon degrees of 
freedom has been used by Berger et al£& to investi- 
gate the Holstcin-Hubbard model at various band fill- 
ings, and Green function QMC simulations for the half- 
filled Holstein model of spinless fermions have been per- 
formed by McKenzie et al£^ Finally, the stochastic se- 
ries expansion MC technique has been applied recently 
to an extended, one-dimensional Hubbard model with an 
electron-phonon interaction of the SSH type^ In con- 
trast to other work, the phonons are treated in second 
quantization. Although the method allows simulations 
on large lattices in one dimension, it relies on an upper 
limit for the number of phonons at each site which makes 
it difficult to study the regime of small phonon frequen- 
cies and/or strong coupling. 

In addition to the specific shortcomings of each method 
such as, e.g., the restriction to a single fermion, or to 
one spatial dimension, or to zero temperature, all pre- 
vious simulations of the Holstein model were limited to 
some extent by autocorrelations. If the phonon degrees 
of freedom are not integrated out analytically, these cor- 
relations predominantly come from the free harmonic- 
oscillator dynamics, especially in the adiabatic regime of 
small phonon frequency. This often leads to an enor- 
mous computational effort even for rather small systems 
and intermediate temperatures. 

In this paper we first present a simple variational 
approach using a generalized form of the Lang-Firsov 
transformation which, in the one-electron case consid- 
ered here, gives surprisingly good results and yields valu- 
able insight into the mechanism of polaron formation. 
The full Hamiltonian resulting from the standard ver- 
sion of the canonical Lang-Firsov transformation is then 
used as the starting point for a QMC method which is 
free of any uncontrolled approximations. Due to the fact 
that the Lang-Firsov transformation contains the crucial 
electronic influences on the phonons, the Monte Carlo 
simulation for the phonon degrees of freedom can be 
based only on the purely phononic part of the trans- 



formed Hamiltonian. The electronic contributions can 
then be allowed for by reweighting of the probability dis- 
tribution, corresponding to an exact treatment of the 
fermion degrees of freedom. This enables us to com- 
pletely ignore the electronic weights in the updating pro- 
cess, and thereby dramatically reduce the computational 
effort. Finally, we introduce a principal component repre- 
sentation of the phonon coordinates, which allows exact 
sampling of the phonons and avoids all autocorrelations. 

The paper is organized as follows. We briefly review 
the Holstein model in Sec.[H] In Sec. IIIII we discuss the 
aforementioned transformations of the Hamiltonian. Sec- 
tion IIVI is devoted to the variational polaron approach, 
while the QMC method for the Lang-Firsov transformed 
model is presented in Sec. Section IVII describes the 
reweighting method, and in Sec. IVIII the representation 
of the phonons in principal components is introduced. 
Results obtained with the given methods are presented 
in Sees. IVIIII and IIXI Finally, Sec. |X] contains our con- 
clusions. 



II. THE HOLSTEIN MODEL 

The Holstein model has been introduced in the 1950's^ 
and describes a system of tight-binding conduction elec- 
trons coupled to a dispersionless phonon mode. If we 
express the phonon operators in terms of their natural 
units, the Hamiltonian takes the form 

H = K + P + I, 
K = ~t c «7 C i CT > 

i 

I = -ayjuxj. (1) 

i 

Here cJ CT (c !(T ) creates (annihilates) an electron of spin 
a at lattice site i, Xi and pi denote the displacement 
and momentum of a harmonic oscillator at site i, and 
hi = J2 a fii<T with hi a — c\ G Ci a . The last term / describes 
the local coupling of the dispersionless Einstein phonons 
to the electron density hi. In the first term, the sym- 
bol (ij) denotes a summation over all nearest-neighbor 
hopping pairs and The parameters of the 

model are the hopping integral i, the phonon energy u>, 
and the electron-phonon coupling constant a. We intro- 
duce the commonly used dimensionless coupling constant 
A = a 2 /(uiW), where W — Atd is the bare bandwidth in 
d dimensions. We also define the dimensionless phonon 
frequency u> = u/t and express all energies in units of 
t. Thus the model depends on two independent parame- 
ters, ui and A. Throughout this paper periodic boundary 
conditions in real space are assumed. 

The Holstein model has been investigated intensively 
in the past, using a large variety of methods. Due 
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to the large amount of literature available we restrict 
the discussion to the case of a single electron in an 
otherwise empty lattice with which this paper is con- 
cerned. The latter is generally known as the Holstein po- 
laron problem and still constitutes a complicated many- 
body problem. Weak-coupling perturbation theory has 
been found to be accurate only for very small coupling 
strengths A when the phonon frequency is low^& In 
the strong-coupling regime, an adiabatic small-polaron 
approximation^^ has been found to work extremely well 
for small values of u> (Ref. l3Sh . while a perturbation the- 
ory based on the Lang-Firsov transformation^ 9 , gives ac- 
curate results for ui ^> 1 (Refs. 1.361 and l38f) . Discrepan- 
cies remain, however, in the regime of intermediate cou- 
pling and phonon frequency4& To bridge this gap, a lot 
of numerical work has been done using exact diagonal- 
ization (ED) methods, density-matrix renormalization- 
group (DMRG) studies, QMC methods and variational 
methods. ED is limited in the accessible parameter 
range, since it requires a truncation of the Hilbert space 
associated with the phonon degrees of freedom. With in- 
creasing electron-phonon coupling strength, for example, 
more and more phonon states have to be included to ob- 
tain converged resultsi 36 i 38 i 41 ' 42 i 43 ' 44 i 45 i 46 i 47 which makes 
it difficult to study clusters of reasonable size in the 
strong or even intermediate coupling regime, especially 
for small phonon frequencies. At this point calculations 
based on DMRG set in, which use an optimized phonon 
basis to reduce the size of the Hilbert space 

48.43.50.51 _A_ n 

other possible approach are the so-called cluster meth- 
ods which exploit exact information on small clusters 
to obtain approximate results for infinite systemsi^S 
Moreover, a number of variational methods have been 
developed which give very accurate results over a wide 
range of parametersi 54 i 55 i 56 i 57 i 58 i 59 i 60 i 61 i 62 As discussed 
in Sec. |]J various QMC methods have been developed 
for the Holstein model. The polaron problem consid- 
ered here has been investigated by Hirsch et a/. ) 13 i 14 De 
Raedt and Lagendijkj 2 ^ 2 ^^ KornilovitchjS^Kornilovitch 
and Piker 8 . Kornilovitchri 3 ^ and Mishchenko et alw^ 
Finally, the Holstein polaron has also been studied in 
the infinite-dimensional limit using dynamical mean-field 
theory** 3 . 

For the one-dimensional case, on which we will focus 
here, the general picture emerging from these investiga- 
tions is as follows (see, e.g., Ref. |5(J). Starting from the 
noninteracting system (A = 0) the electron is gradually 
dressed with a coherent multi-phonon cloud as the cou- 
pling increases. For A < 1 and X/u) < 0.5 the result- 
ing quasiparticle remains in a Bloch-like state, with the 
phonon cloud giving rise to an increased effective mass. 
In the strong-coupling regime (A > 1 and A/a) > 0.5) 
the electron becomes self-trapped by the induced lat- 
tice distortion and this object — trapped electron plus 
distortion — is usually called a small polaron. The tran- 
sition from weak to strong coupling is continuous^ and 
the term "large polaron" is often used to describe an elec- 
tron dressed with a phonon cloud extending over more 



than one lattice site. The polaronic effects become more 
dominant as the phonon frequency approaches the adia- 
batic limit Q — > 0. In QMC simulations, small values of u> 
introduce two very different time scales for the evolution 
of electrons and phonons, respectively. This gives rise to 
the problems mentioned above and, in fact, many QMC 
simulations have been done only for Q > 1. 



III. EXTENDED LANG-FIRSOV 
TRANSFORMATION 

The canonical Lang-Firsov transformation 3 ^ has been 
used extensively to study the polaron problem. A well- 
known approximation due to Holstein^ consists of re- 
placing the transformed hopping term by its expectation 
value with respect to a zero-phonon state, thus neglect- 
ing phonon emission and absorption during the hopping 
process. This approximation, which we shall call the 
Holstein-Lang-Firsov (HLF) approximation, was found 
to give reliable results only in the strong-coupling and/or 
nonadiabatic limit A,ui > 1 fRefs. IHol liiil and Ej) . More 
refined approaches based on strong-coupling perturba- 
tion theory provide an accurate description of the Hol- 
stein polaron over a large range of parametersi 3 ^ 3 ^ In the 
limit A = oo, the hopping term in Hamiltonian £Q can be 
neglected, and the Lang-Firsov transformation allows an 
exact solution of the resulting single-site problem. 40 The 
transformation has also been used in combination with 
numerical methods^Si 6 ^ However we are not aware of 
any QMC simulation based on the transformed model. 

We define the unitary operator 

U = e v , v = i ^2 lijfiiPj , ( 2 ) 
ij 

where i and j run over lattice sites, and with real pa- 
rameters jij. U as defined in Eq. @ has the form of a 
translation operator. Given an electron at lattice site i, it 
corresponds to a displacement of the harmonic oscillators 
at all sites j, j = 1, . . . ,N, by 7^ . Hence the transfor- 
mation describes a nonlocal phonon cloud surrounding 
an electron. This corresponds to the well-known concept 
of a large polaron, which extends over more than one lat- 
tice site. Using the transformation O = U O U* we find 
for the transformed operators 

Xi Xi -\- ^ ^fijUj 5 Pi Pi 

3 

4 = cle^i^ , c^c^e" 1 ^^ . (3) 
Inserting these results into Eq. the transformed 
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Hamiltonian becomes 

H = K + P + I ep + I ee 

I ep = E njXj(wyij - aSjj) 




Here the term I ep describes the coupling between elec- 
trons and phonons, while I cc represents an effective 
electron-electron interaction. The Hamiltonian @ will 
be the starting point for the variational polaron approach 
presented in the following section. 

A more suitable approach for QMC simulations, how- 
ever, is given by requiring that the electron-phonon terms 
cancel. Then 7^ = 7^ with 7 = yj XW/uj and we ob- 
tain the standard Lang-Firsov transformation with the 
transformation operator 

U a = e u ° , u = i7 ^2 hipi . (5) 

i 

In contrast to the extended polaron cloud, defined by 
Eq. J2J ! now only the oscillator at the site of the electron 
is affected. The transformed Hamiltonian reads 

H = K + P + Q, 

k = -* E A^J" 1 ^-^ , 

Q = -i 7 2 -E^- (8) 

i 

In the HLF or small-polaron approximation, the ground 
state of the transformed Hamiltonian is approximated 
by leaving all phonons in the ground state. It has been 
shown^i that the small-polaron wave function becomes 
exact in the strong-coupling, nonadiabatic limit, and 
agrees qualitatively with the exact results also in the 
intermediate coupling regime. As discussed by Zhang 
et alr^ the HLF approximation gives an overestimated 
shift of the equilibrium position of the oscillator in the 
presence of an electron, and does not reproduce the re- 
tardation effects when an electron hops onto a previously 
unoccupied site. Nevertheless, the local lattice distortion 
at the site of the electron contains the crucial impact 
of the electron on the lattice. Consequently, the trans- 
formed Hamiltonian © should be a good starting point 
to perform QMC simulations, which merely need to sim- 
ulate small fluctuations around the zero-point motion. 
Indeed we will see in Sec. IIXI that the expectation val- 
ues of the phonon operators in the transformed Holstcin 
model remain close to the results of the free-oscillator 
case over the whole range of the electron-phonon cou- 
pling. This makes sampling of the phonon degrees of 



freedom much more efficient. In principle, it would also 
be possible to develop a QMC algorithm starting with 
Hamiltonian Q, with the parameters 7^ being deter- 
mined by the variational method discussed in the follow- 
ing section. However, we will see that the simple (local) 
Lang-Firsov transformation is already sufficient to obtain 
a very efficient QMC method. 

From Eq. J§} it is obvious that the standard Lang- 
Firsov transformation on the one hand removes the 
electron-phonon coupling term, but on the other hand 
introduces complex valued hopping integrals which de- 
pend on the phonon momenta at the lattice sites in- 
volved in the hopping process. Moreover, for more than 
one electron in the system, the last term Q introduces a 
Hubbard-like attractive interaction. In the case of the ex- 
tended transformation, the electron-phonon interaction 
term cannot be eliminated entirely, the hopping term in- 
volves all phonon momenta pi as well as the parameters 
7y , and the electron-electron interaction becomes long 
ranged. For these reasons it is expedient to base the 
QMC simulation on the local Lang-Firsov transforma- 
tion ©. 

For simplicity, we restrict ourselves in the present 
study to the case of a single electron so that hihj = hiSij, 
The electron-electron interaction term in Hamiltonian lf4"| 
becomes 

lee = E " l E 7 ^ ~ "^M ' ( 7 ) 

while the corresponding term in the local Lang-Firsov 
transformation [last term of Hamiltonian ©] reduces to 

Q -XW/2 = -E P . (8) 

Equations Q and JHJ) both describe a shift in energy 
resulting from the original electron-electron interaction 
which is usually called the polaron binding energy Ep. 

IV. VARIATIONAL POLARON APPROACH 

Here we present a simple variational method which is 
based on the extended transformation discussed in the 
preceding section. Similar work along these lines using 
different transformations of the Hamiltonian as well as 
physically motivated wave f unction s can be found, for 
example, in Refs. IHsl and l6fl6]l As discussed above, 
the zero-phonon ansatz of the simple HLF approxima- 
tion gives reliable results only in the limit of large ui and 
A. Whereas in HLF the parameter 7 of the Lang-Firsov 
transformation is chosen such that the electron-phonon 
coupling term J ep vanishes, in the variational polaron ap- 
proach ( VPA) , we treat the 7^ as variational parameters 
which are determined by minimizing the ground-state 
energy in a zero-phonon basis. Like the HLF approx- 
imation, the VPA becomes exact in the weak-coupling 
limit A — > and in the nonadiabatic strong-coupling limit 
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A,w — > oo. We will see in Sec. IVIIII that the VPA also 
gives very accurate results for large phonon frequencies, 
Q 1. This can easily be understood keeping in mind 
the discussion of the validity of the HLF approximation 
given in the preceding section. While the HLF ansatz 
overestimates the displacement of the local oscillator in 
the presence of an electron, the VPA determines this shift 
variationally. Moreover, the missing retardation effects in 
the response of the oscillator to an electron hopping onto 
the site become negligible for large phonon frequencies. 
Therefore, in addition to the cases stated above, the VPA 
also becomes exact in the nonadiabatic limit Co — ► oo. 
Although the limitations of the VPA in or near the adia- 
batic regime will clearly emerge when we discuss results 
in Sec. lVIIll it works surprisingly well if we keep in mind 
the simplicity of the method. Moreover, the reasons for 
the failure of the VPA in certain parameter regimes are 
physically clear and can easily be interpreted. 

For translationally invariant systems the displacement 
fields satisfy the condition 7^ = J\i-j\- Inserting this 
relation into Eq. (7J, the expression inside the brackets 
becomes independent of the index i. For the single elec- 
tron case with hi = 1 we have 



!cc = \ E 
1 



«7o 



We solve the eigenvalue problem of the transformed 
Hamiltonian in a zero-phonon basis for which we 
make the ansatz 

IO=cLlO)®ni^ ) ) I 1 = 1,- -.,N, (10) 



where \<fi^} denotes the ground state of the harmonic 
oscillator at site v. For simplicity, we restrict ourselves 
to one dimension, although the method can easily be ex- 
tended to higher dimensions. The matrix elements of the 
hopping term in this basis are 



(l\K\l') 



= -tw Yl / & x 4>{% + llu)^{x + 7i/„) 

= -t w e-\^^-i v+l -v)\ (11) 

where tw = td^u^ is nonzero for nearest-neighbor hop- 
ping pairs I' = I ± 1 and 4>(x) is the harmonic-oscillator 
ground-state wave function in coordinate space. The ma- 
trix elements of the other terms of Hamiltonian are 
easily evaluated yielding 



(i\p\n 
(i\i cp \n 
(i\i cc \n 



= Sw 



0, 



Sw 



2 ' 



E^ 2 



«7o 



(12) 



The eigenstates of the transformed Hamiltonian Q in 
the zero-phonon subspace, spanned by the basis states 
defined in Eq. Hl()[) . are 



iik>=4jo>®rw 



(13) 



with energy 
E(k) 



£k + ? V+-^r 7; 2 -a7o 



-t J2 e ikS e-^ 4 ^,A^~f»+s)\ (14) 
s=±i 



where E k denotes the kinetic energy of the electron. 
Defining the Fourier-transformed parameters "f q as 



and using (7; € M.' 



E %^i eiqS = E % 2 cos qS 

q q 



(15) 



(16) 



(9) the kinetic energy can be written as 



£ k = _ t y^ e ife5 e -(l/2)E a (l-cos g ,5)7| 
5 

= e (k)e-^ 2 ^^- cos ^ 
= £off(fc) , 



(17) 



where ?o(fc) = — 2<cosfc is the tight-binding dispersion 
in one dimension. Using these results the ground-state 
energy finally becomes 



I 2 

q 



Ev ( 18 ) 



The variational parameters j p are determined by 

dE a \ 

TTz- = -7 P e c ff(p)(l - cosp) + u% 7=F?=°- ( 19 ) 

C7 P y/N 



The values for 7 P which minimize the energy E can then 
be obtained from 



1 



7 P 



N u> + e c ff(p)(l - cosp) 



(20) 



As e e ft depends on the set of parameters 7 P , this equation 
has to be solved self-consistently. Equation l|20|l has a 
typical random-phase approximation form, which is rea- 
sonable since a variational ansatz for the wave function 
of the untransformed Hamiltonian can be written as [see 
also Eq. ©] 



Uj = 7f E ^ c l e- 1 ^™* |0> ® II 14' 



(21) 
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In addition to the total energy given by Eq. I|18(l . we are 
also interested in the quasi-particle weight for momentum 
k = 0, defined as 



zq = (0| c k=Q:<T \4>o) 



(22) 

Here |i/'o) denotes the ground state with one electron of 
momentum p = and the oscillators in the ground state 
|0o). Fourier transformation leads to 



l^(0o|(O|^c] CT |O)|0 o ) 
^(0o|e- iS ^ Wfc |0o) 



(23) 

where we have used the same steps as in Eq. <|11|) . 

Results obtained with the VPA will be presented in 

Sec. lvTnl 



V. MONTE CARLO FOR THE TRANSFORMED 
MODEL 

In contrast to the approximate variational approach 
presented in the preceding section, the QMC method 
discussed here is based on the exact Lang-Firsov trans- 
formation of the Holstein Hamiltonian. Therefore, the 
method is exact apart from statistical errors and Trotter 
discretization, as discussed in Sec. H] 



A. Partition function 

We begin with the evaluation of the partition function 
Z = Tre~ 0H = Tver 1311 , with H given by Eq. ©. As 
indicated in the preceding section, for the case of a single 
electron, the last term in Hamiltonian © represents a 
constant energy shift. Moreover we can drop spin indices 
and are left with the Hamiltonian 

H = K + P-Ep. (24) 

The polaron binding energy given by Eq. JHJ can be ne- 
glected during the QMC simulation, and needs only to be 
considered in calculating the total energy. For simplicity, 
we only consider the one-dimensional case here, although 
the generalization to higher dimensions is a simple mat- 
ter. Using the Suzuki- Trotter decomposition we obtain^ 



-f3H 



(e -ArK e -ArP Pe -ArP^L 



(25) 



where (3 = (ksT)- 1 and At = 0/L. Splitting up the 
trace into a bosonic and a fermionic part and inserting L 
complete sets of momentum eigenstates^ we derive the 
approximation for the partition function 



Z L = Tr 



dpidp 2 • • • dp L (pi | U |p 2 ) • • • (PL | U \pi 



(26) 



where dp T = Yii dpi. T . Each matrix element can be eval- 
uated by inserting a complete set of phonon coordinate 
eigenstates J dx T \x T )(x T \. All x T integrals are of Gaus- 
sian form and can easily be carried out. The result is 



C = 



2tt 
ujAt 



(27) 



The normalization factor in front of the exponential has 
to be taken into account in the calculation of the total 
energy but cancels when we measure other observables. 
With the abbreviation T>p = dpidp 2 ■ ■ • dp^ the partition 
function finally becomes 



Z L = C 



NL 



Vp w\, Wi 



(28) 



with the abbreviations 

w h = e- ATS \ io f = Tr f fi, Q = Y[ e~ Ari?T . (29) 



Here K T is K with the phonon operators pt replaced by 
the momentum p^ T on the rth Trotter slice. The expo- 
nential of the hopping term may in the single-electron 
case be written as 



e -ArA- T = DtK ^ 



(30) 



Kjf — (e j , [D T )jji — 5jjie 1P3)T , 

v ' 33' 

where h th is the N x N tight-binding hopping matrix. 
Thus we have the same matrix k for every time slice, 
which is transformed by the diagonal unitary matri- 
ces D T . In our one-electron case, the fermionic weight 
Wf = J2 n ( n \ ^ l n ) ^ S gi yen by the sum over the diagonal 
elements of the matrix representation of f2 in the basis of 
one-electron states 



\n) = C t |0) . 



(31) 



The bosonic action in Eq. I|29|l contains only classical 
variables and takes the form 



\ ^ 2 



Yl ( Pi > T ~ P^+i) 2 . (32) 



Sb- 1 2a,(Ar) 

where the indices i = 1, . . . , N and r — 1 , . . . , L run over 
all lattice sites and Trotter times, respectively, with the 
periodic boundary conditions Pi t L+i — Pi.i- It may also 
be written as 



(33) 



with pi — (pi.ii • • • iPi,L) and a "periodic" tridiagonal 
L x L matrix A with nonzero elements 



oj 1 
2 + 



1 

ujAt 2 



(34) 
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Since 2x is a trace, it follows that A\ l = = 
-1/(uAt 2 ). 

At this stage, with the above result for the partition 
function, a QMC simulation of the transformed Holstcin 
model would proceed as follows. In each MC step, a pair 
of indices (io, tq) on the N x L lattice of phonon mo- 
menta Pi iT is chosen at random. At this site, a change 
Pi ,T l— * Pio-To + <-^P °f the phonon configuration is pro- 
posed. To decide upon the acceptance of the new config- 
uration using the Metropolis algorithm, the correspond- 
ing weights WbWf and w' h w' f have to be calculated. Due 
to the local updating process, the change of the bosonic 
weight Awh — w' h — Wh can easily be obtained. In con- 
trast, the fermionic weight requires the evaluation of the 
L fold matrix product appearing in the definition of SI 
in Eq. 129|) . The numerical effort for the calculation of 
W{ may be reduced by varying To sequentially from 1 to 
L instead of picking random values. In this case the cal- 
culation of the new fermionic weight, after the change of 
a single phonon momentum, can be reduced to only two 
matrix multiplications. Similar to other MC methods, a 
warm-up phase at the beginning of the simulation would 
be required for each set of parameters. An additional 
difficulty arises from the fact that, for the transformed 
model, the fermionic weight Wf is no longer strictly pos- 
itive, even for the case of a single electron in one di- 
mension. This is a consequence of the complex-valued 
hopping integrals, in contrast to simulations of, e.g., the 
Hubbard model, where a minus-sign problem occurs as 
a consequence of the Fermi statistics of the electrons^ 
Here the average sign of Wf is smallest in the regime of 
small phonon frequency and low temperature. The sign 
problem is most pronounced for intermediate values of 
the electron-phonon coupling strength A, where the cross 
over from a large to a small polaron occurs. However, 
in one dimension, it is not severe and reduces with in- 
creasing system size. For example, calculations in the 
most critical regime fit = 10, u) = 0.1, and A ss 1 have 
shown that (sign) = (w{)/(\wf |) increases from about 0.5 
for N = 4 to about 0.85 for N = 16. Nevertheless, it re- 
mains to be seen to what extent the number of electrons 
and the dimensionality of the system affect the situation. 

A related QMC approach to the original Holstein 
Hamiltonian involves a very similar derivation^ to 
obtain the partition function, also in the one-electron 
limit. In fact the bosonic action St, takes exactly the 
same form, with p replaced by x. The main difference 
is the fermionic part of the partition function, contained 
in the matrix Q. While the Lang-Firsov transformation 
leads to a complicated hopping term, the standard ap- 
proach for the untransformed model only includes the 
bare hopping operator given by Eq. Q). However, an in- 
teraction term / [cf. Eq. appears, which contains the 
phonon coordinate x. Hence the matrix Jl is replaced by 

L 

V = \{kV t , (Vr) n , = S n , e A — (35) 

T=l 



and the path integral in the partition function [Eq. 128|) ] 
is over all coordinates x instead of the momenta p. Apart 
from the fact that the coordinates x are sampled instead 
of the phonon momenta, the QMC procedure for the 
untransformed model is identical to the simulation de- 
scribed above. We shall refer to this less sophisticated 
QMC method for the original Holstein Hamiltonian as 
the standard approach. For A = a = 0, i.e., no electron- 
phonon coupling, we have a set of N independent har- 
monic oscillators, and both approaches are alike. 



B. Problems with the standard approach 

Let us briefly consider the noninteracting limit, in 
which the partition function can be written as Zl ~ 
fT>pe~ ATSb . As discussed by Batrouni and Scalettar^ 
the difficulties encountered in QMC simulations, even 
for the simple case of a single (N — 1) harmonic os- 
cillator, arise from the large condition number, i.e., the 
ratio of largest to smallest eigenvalue, of the bosonic ac- 
tion Sb- For small values of At this ratio is propor- 
tional to {loAt)~ 2 (Ref. I66I) . leading to autocorrelation 
times which grow quadratically with decreasing phonon 
frequency and the number L of Trotter times. The phys- 
ical reason for these correlations becomes obvious if we 
look at the bosonic action [Eq. I|32l) ]. The latter can be 
thought of as being proportional to the energy of a given 
phonon configuration, E = ArSb- While the first term 
corresponds to the kinetic energy of the oscillators, the 
second term describes a coupling in imaginary time, i.e., 
a pure quantum effect. As pointed out by Batrouni and 
Scalettar^ large changes of a single phonon degree of 
freedom, pi tT say, are very unlikely to be accepted due to 
the energy change proportional to l/(wAr), which arises 
from the coupling to Pi, T ±i- However, a QMC simula- 
tion with only small local changes is extremely ineffective 
in sampling the relevant regions of phase space. There- 
fore, successive phonon configurations will be highly cor- 
related. A possible solution might be the use of global 
updating schemes. Alternatively, the situation could be 
improved by transforming to the normal modes of the 
phonons, so that different step sizes can be used in up- 
dating each mode. We will see in Sec. IVIII that the prin- 
cipal component representation can indeed be used to 
completely eliminate these difficulties. 

In addition to the above-mentioned autocorrelations, 
which are in fact independent of any electronic influ- 
ences, standard simulations of the Holstein model be- 
come very difficult in the regime where polaron effects are 
large. This occurs at low temperatures, small phonon fre- 
quencies Co < 1, and for intermediate or strong electron- 
phonon coupling A > 1. Unfortunately, these are ex- 
actly the parameters of interest for simulations of real 
substances such as, e.g., the manganitesi^ To illustrate 
the physical origin of these problems let us consider the 
case of a single electron in the Holstein model. As dis- 
cussed in Sec. [Hj in the polaronic regime, the electron 
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drags with it a cloud of phonons which corresponds to a 
more or less localized lattice distortion. When the elec- 
tron hops from site A (with a displaced oscillator cor- 
responding to a small polaron), say, to a neighboring, 
previously unoccupied site B (with the oscillator in its 
undisplaced ground state) during a QMC simulation, the 
current phonon configuration is no longer energetically 
favorable. Clearly, the oscillator at site A has to return 
to its undisplaced ground state, while a corresponding 
phonon cloud has to be built up at site B. Such distor- 
tions of the lattice in the presence of an electron are large 
compared to the zero-point motion of the oscillator. On 
the other hand, only small changes of the current config- 
uration will be accepted in the simulation. Consequently 
it takes an enormous number of single updates to ob- 
tain the new configuration in which the polaron has com- 
pletely moved to site B. Obviously these polaron effects 
also give rise to strongly autocorrelated configurations, 
thereby dramatically increasing the numerical effort for 
the simulation. These problems due to polaron formation 
can be overcome by using the Lang-Firsov transformed 
model. The transformation separates the large displace- 
ments of the local oscillators, due to polaron effects, from 
the free-oscillator dynamics which correspond to vibra- 
tions around the shifted equilibrium positions. The quan- 
tities to be sampled, namely the phonon momenta p, only 
show a weak dependence on the electron-phonon coupling 
strength A, in stark contrast to the coordinates x in the 
original, untransformed model, whose expectation val- 
ues grow linearly with A in the strong-coupling regime. 
In fact, the QMC results obtained for the transformed 
model (see also Sec. IIXII show that the statistical er- 
rors increase in the intermediate coupling regime A w 1, 
but decrease again as we approach the strong-coupling 
limit. This is in perfect agreement with the fact that the 
the Lang-Firsov transformation diagonalizes the Hamil- 
tonian in the strong-coupling or atomic limit A — > oo 
(see Sec. IHIfl . so that the QMC method based on the 
transformed model becomes more and more efficient as A 



energy which is defined as 



increases. 



Observables 



Thermodynamic expectation values 
(O) = Z- 1 Tr 6 e~ pH = Z' 1 Tr 6 



(36) 



of observables O are computed in the Lang-Firsov trans- 
formed representation via 

(O) = Z- 1 Tr f / dp (p\ 6 e-W \p) . (37) 



In this paper we are interested in the kinetic energy of the 
electron, the total energy, the mean square of the phonon 
momenta, and the momentum distribution n(k) = (ctefe) 
for various wave vectors k. We begin with the kinetic 



E k = (K) = -tZ- 1 Tr ( 4 c i e'ft*-^ e^ H ) . 

Using the same steps as in the derivation of the parti- 
tion function (see Sec. IV A|l . and absorbing the additional 
phase factor in a matrix M — d\VLD\ [see Eq. <|3(J|) ]. we 
find 

E k = -tZl x Y^ I Vpw h ^2{n\Mc\c 3 \n) 

m 

= -tZl 1 I Vpw h (j\M\t) 

with one-electron states \n) as defined in Eq. (|31|l . Using 
the matrix elements My = (i\M \j) and the expectation 
values 



<0) b = 



/ Vpw h 0(p) 



(39) 



/ Vpw h 

with respect to the purely phononic weights Wh we obtain 



E k = -t 



E< (Mi 



ii/b 



(40) 



Here we have already taken into account the reweighting 
method which will be discussed in detail in the following 
section. The total energy can be obtained from the ther- 
modynamic relation E = —d(hiZ)/df3, with Z given by 
Eq. QgJ. The result is 



(41) 



EL 



P h 



N 



2At 2wAt 2 L 



where Ep is defined in Eq. © and the expectation values 
are calculated according to Eq. (|4^|) given below. To 
compare with other work we subtract the ground-state 
energy of the phonons, So,ph = Nui/2. Finally, n(k) can 
be obtained using Fourier transformation. The result is 



n{k) 



1 E» (M; 



ij V"*7/b' 



,ik(i-j) 



N 



T,i(M; 



ii/b 



(42) 



with k from the first Brillouin zone and the same matrix 
M as in the case of the kinetic energy. 



VI. REWEIGHTING 

In typical QMC simulations a large amount of the to- 
tal computational effort goes into the calculation of the 
probability for the acceptance of a proposed change of 
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the configuration. This probability is usually determined 
by the ratio of the weights of the new and the old con- 
figuration, as in the Metropolis algorithm used here. In 
the notation of Sec. [V] this involves the calculation of 
u»b and ui{ for the two configurations, S and S' say, in 
every MC step. While the change in the bosonic weight, 
Wb(S')/wb(S), is easily calculated for the case of local up- 
dating, the fermionic weight given by Eq. (|29|) involves an 
L-fold matrix product of N x N matrices for each config- 
uration. Although the numerical effort of the evaluation 
of such a matrix product can be reduced by scanning se- 
quentially through the time slices (see Sec. IV A|) it still 
requires a lot of total computer time. 

This can be avoided by reweighting of the probability 
distribution to be sampled. In the case considered here, 
this corresponds to taking into account only the change 
wt,(S')/wb{S) in the bosonic weight, and compensating 
for this by dividing the resulting expectation value by the 
expectation value of the fermionic weight Wf , as has been 
used already in Eq. I4UI leading generally to ratios of the 
form 



(O) = 



(Owfh 

(tOf)b 



(43) 



where the subscript "b," defined in Eq. I|39|) . indicates 
that the average is computed based on Wb only. Fol- 
lowing this procedure, the fermionic weight is treated as 
part of the observables. The splitting into weight Wb and 
observable Ow{ is sensible as long as the variance of Wf 
and Ouif is small, which is the case after the Lang-Firsov 
transformation. This approach has several additional ad- 
vantages. With the reweighting method, the updating of 
the system does no longer require the calculation of Wf in 
every step, but only when measurements are performed. 
Compared to the usual QMC procedure described in 
Sec. IV Al this can save an enormous amount of computer 
time, allowing such simulations to be run on a stan- 
dard personal computer instead of a high-performance 
supercomputer. Additionally, since the updating does 
no longer involve any electronic contributions, it becomes 
independent of the electron-phonon coupling strength A. 
This allows the simultaneous measurement of observables 
for a whole set of values of A in a single MC run. For a 
given phonon configuration, the fermionic weight and the 
observables are measured and stored for each value of the 
coupling. This procedure is repeated until the required 
number of measurements has been made. At the end of 
the simulation an appropriate analysis of the measured 
values is made independently for each A. In contrast, 
the QMC procedure without reweighting (see Sec. IV A|) 
would require a separate run for each value of A, including 
a warm-up phase to equilibrate the system for the current 
set of parameters. We will see in Sec. IVIll that in combi- 
nation with the principal component representation, the 
phonon momenta p can be sampled exactly, removing 
all autocorrelations. This avoids a warm-up phase, and 
measurements can be made after every Monte Carlo step. 
In this final, very efficient procedure, the calculation of 



10, co = 0.1 

5,(0 = 0.1 
5,(0 = 2.0 
5, (0 = 2.0, w/o LF 
5,(0 = 0.1 




FIG. 1: Kullback-Leibler number /ikl as a function of 
electron-phonon coupling A for various sets of the parame- 
ters N (number of sites), f3 (inverse temperature) and Q. As 
indicated, the results for the untransformed model, denoted 
in the legend as "w/o LF," have been scaled by a factor 0.25 
(see text). Error bars are smaller than the symbols shown, 
and lines are guides to the eye only. 



Wf for measurements remains, and is then the most time- 
consuming part of the calculation. Finally, we want to 
point out that, with the use of the reweighting method, 
the electronic degrees of freedom are treated exactly, i.e., 
they are not sampled in the course of the simulation. 

As mentioned in Sec. IV Al the weight Wf for the trans- 
formed model is no longer strictly positive, so that it 
cannot be interpreted as a probability. The usual way to 
deal with such a sign problem is to split the weight into 
Wf = | Wf | sgnwf. Then |wf| can be used as the weight of 
a given configuration in the updating process, while the 
sign is absorbed in the observables. The difference to the 
reweighting method presented here is that instead of the 
sign of Wf, we treat the whole weight Wf as part of the 
observables. 

Despite the obvious advantages of this approach, it 
is necessary to scrutinize whether reweighting does not 
lead to prohibitive statistical noise. If, for example, there 
was too small an overlap of the actual probability distri- 
bution with the one we are sampling with, the method 
would fail. In fact, our calculations have shown that for 
the imtransformed model the reweighting method would 
fail at low temperatures and for critical values of the pa- 
rameters ui and A. 

The distance between two arbitrary probability dis- 
tributions 4>i(x) and foix), each depending on a set of 
variables x, can be measured by the so-called Kullback- 
Leibler number /xkl which is defined as& 



dx 4>i(x) In 



<t>i{x) 
<fa{x) 



(44) 



For 



t>x = cf>2 we have /^kl = 0, while for (f>i ^ c/> 2 
/ikl > 0. The fact that /zkl is a reasonable measure 
for the distance of two distributions is best illustrated 
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by considering two Gaussian deviates <px, 4>2 with vari- 
ance cr 2 , centered at X\ and X2, respectively. In this case 
A'kl = [x% — X2) 2 / {2a 2 ) . For \x%— x%\ = \/2<r, where the 
two peaks begin to be distinguishable, we have ^kl = 1, 
while a large value of //kl — 10, for example, corresponds 
to well-separated Gaussian distributions. Here we use the 
Kullback-Leibler number to investigate the applicability 
of the reweighting method. As long as the Kullback- 
Leibler number is less than or comparable to 1, reweight- 
ing works well, while a Kullback-Leibler number strongly 
exceeding unity indicates severe problems. Two relevant 
distributions in our case are given by 4>i{p) = Wb(p)/Zb 
and 4>2(p) = Wb(p)\w[(p)\/ Zbf, depending on the phonon 
configuration p (or x in the case of the untransformed 
Holstein model). Z^ and Z\>{ are the normalization fac- 
tors of the probability densities <fii(p) and 4>2 (p), and 
Wf has been replaced by its absolute value due to the 
aforementioned sign problem. Inserting these definitions 
into Eq. we find /i K l = hi(|wf|)b — (In |iUf|)t,. Fig- 
ure shows results for /zkl for different parameters j3, 
ui, and N. For A = 0, Wf is independent of the phonon 
configuration so that /ukl = 0. With increasing electron- 
phonon coupling, the difference between the two distri- 
butions becomes larger. For an intermediate value of the 
electron-phonon coupling strength, /ukt takes on a max- 
imum and approaches zero again in the strong-coupling 
limit A — > 00. This is exactly the behavior we would ex- 
pect for the Lang-Firsov transformed model. For A = 
the transformation has no effect and Wf is a constant, just 
as in the case of the untransformed model. In the inter- 
mediate coupling regime, the small-polaron picture me- 
diated by the transformation is not correct as we have an 
extended (large) polaron in this region. However, as the 
coupling increases further, the polaron becomes smaller 
and for A = 00 it is known that the Lang-Firsov transfor- 
mation diagonalizes the Holstein Hamiltonian QJ. The 
dependence of /j,kl on the temperature and the phonon 
frequency is also in perfect agreement with the physi- 
cal picture of the Holstein polaron. As fit increases, po- 
laron effects become more prominent. The same effect 
occurs if we decrease Co, and in both cases the maximum 
of hkl increases. In Fig. ^ the result for a system of 
eight lattice sites is also shown. The maximum in /^kl 
is clearly smaller than for the four-site cluster. Calcu- 
lations for even larger clusters (not shown) reveal that 
the maximum in /l*kl decreases further indicating that 
the overlap between <f>\ and <f>2 increases as N — > 00. 
This behavior agrees well with the influence of finite- 
size effects in the transition region as pointed out before 
by Marsiglio^ As the system size increases, the cross 
over becomes smoother in agreement with the fact that 
the ground state of the Holstein polaron is an analytic 
function of the coupling A (Ref. 1641) . This point will be 
further illustrated in Sec. IIXI To summarize, for all pa- 
rameters shown in Fig.^ the maximum of /ljkl lies below 
/iKL ~ 1, so that we can conclude that the two distribu- 
tions are indeed very close and the reweighting method 
can be successfully applied. 



We have also calculated the Kullback-Leibler number 
for the case of the untransformed model, denoted in Fig.Q] 
as "w/o LF," for which |ujf| = Wf. The result has been 
divided by a factor 4 to allow a better representation 
in Fig. The difference between <j>\ and 4>2 increases 
strongly with A and reaches large values of /ikt > 10 
already in the intermediate coupling regime 1 < A < 
2. Hence we cannot expect the reweighting method to 
work in this case. Finally we want to point out that the 
distance between <pi and 4>2 may not affect all observables 
in the same way. A detailed analysis for each observable 
O would be based on the Kullback-Leibler distance of the 
marginal probability densities 

p a (o) = J dx p(o\x) p a (x) = J dx 5(o - 0(x)) p a (x) , 

where 0(x) is the value of the observable for a given 
configuration x and a = 1,2 for the two distributions 
under consideration. 

In summary, the reweighting method, together with 
the Lang-Firsov transformation, allows us to sample a 
system of independent oscillators, while all the influence 
of the electrons is transferred to the observable, thereby 
strongly reducing the numerical effort. In order to ob- 
tain a reliable error analysis for observables calculated 
according to Eq. (|4*31) . the jackknife procedure^ has been 
applied. 

VII. PRINCIPAL COMPONENT 
REPRESENTATION 

Although the reweighting method allows us, in princi- 
ple, to skip enough sweeps between measurements to re- 
duce autocorrelations to a minimum, the computational 
effort for these Monte Carlo updates can become the most 
time-consuming part of the simulation. Even though a 
single phonon update requires negligible computer time 
compared to the evaluation of the fcrmionic weight, in 
the critical parameter regime, an enormous number of 
such steps will be necessary between successive measure- 
ments. Moreover, reliable results can only be obtained 
when long enough Monte Carlo runs have been performed 
to see even the longest autocorrelation times. In this 
section, we present a principal component representation 
for the phonon degrees of freedom, which enables us to 
create completely uncorrelated samples of phonon con- 
figurations. 

In order to illustrate the severe problem of autocorre- 
lations with standard updates of the phonons, we have 
calculated the integrated autocorrelation time rf nt for the 
phonon momenta p. r; nt is a direct measure for the num- 
ber of MC steps which have to be skipped between mea- 
surements in order to obtain uncorrelated data, and is 
usually given in units of sweeps. We define a sweep as N 
times L proposed local changes of the phonon configura- 
tion. For a four-site system, for example, with fit = 5, 
A = 2, Q = 2, and At = 0.05 we find tL ~ 500. This 
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corresponds to an autocorrelation time of about 2 x 10 5 
single MC steps. For smaller phonon frequencies, r- mt 
increases strongly. For Q = 1 and the same At, the au- 
tocorrelation time is already ~ 1700 sweeps, which agrees 
quite well with the (uAt) -2 dependence of the correla- 
tions for A = given in Sec. IV Bl The dependence of 
r kit on t ne coupling strength A is relatively weak, and we 
have found no systematic behavior of T? t as a function 
of A. Depending on the other parameters, the autocor- 
relation times were observed to increase or even decrease 
slightly as A is increased. This behavior can be ascribed 
to the exact treatment of the fcrmion degrees of free- 
dom. As we are not sampling the hopping process of the 
single electron considered here, no autocorrelations due 
to the resulting reaction of the harmonic oscillators to 
the electronic motion (see Sec. IVBfl can occur. More- 
over, even if we would sample the electronic degrees of 
freedom in the QMC simulation, these autocorrelations 
would still be strongly reduced as long as we use the 
Lang-Firsov transformed model. This is a consequence 
of the fact that the large displacements of the oscillators 
in the presence of an electron are explicitly contained 
in the Hamiltonian © . Finally, as the number of lattice 
sites is varied, T? nt remains constant in units of sweeps for 
our single-electron simulations. We also determined the 
autocorrelation times for observables such as, e.g., the 
kinetic energy. Although Tint is smaller for electronic ob- 
servables, the problem still exists, and the determination 
of the autocorrelation times for various parameter sets is 
vital to obtain reliable results. This usually requires very 
long QMC runs and a lot of CPU time. 

As indicated in Sec. IV Bl the autocorrelations which 
arise from the structure of the bosonic action St, [see 
Eq. I|32|)] may be overcome by a transformation to the 
normal modes of the system. Here we represent the 
bosonic action St, in terms of its normal modes along 
the imaginary time axis. This allows us to sample com- 
pletely uncorrelated phonon configurations. In combina- 
tion with the reweighting method the fermion degrees of 
freedom are treated exactly, so that our QMC method is 
indeed free of any autocorrelations. This greatly simpli- 
fies calculations, since it makes the usual binning analysis 
(to determine the autocorrelation times) obsolete and, 
more importantly, leads to significantly shorter simula- 
tion times. 

All this can be achieved with the simple but effective 
idea of a transformation to principal components (PCs). 
To this end let us recall the form of the bosonic action 
given by Eq. i|33|) which can also be written as 

S b = 5> T 4* = £p^ V2 ^ /2 Pi :£«, T 'G (45) 

i i i 

with the PCs £j = A 1 / 2 pi 1 in terms of which the bosonic 
weight takes the simple Gaussian form 

w h = e ~ AT ^ (46) 

The QMC can now be performed directly in terms of 
the new variables £. To calculate observables we have 



to transform back to the phonon momenta p using the 
matrix A^ 1 / 2 . Comparison with Eq. H33[) shows that in- 
stead of the ill-conditioned matrix A we now have the 
ideal structure that we can easily generate exact samples 
of a Gaussian distribution. In terms of the new coor- 
dinates £, the probability distribution can be sampled 
exactly, e.g., by the Box-Muller method^ In contrast to 
a standard Markov chain MC simulation, every new con- 
figuration is accepted, and measurements of observables 
can be made at each step. 

From the definition of the PCs it is obvious that an 
update of a single variable £j )T) say, actually corresponds 
to a change of all phonon coordinates pi yT i , t' = 1, . . . , L. 
Thus, in terms of the original phonon coordinates pi, the 
updating loses its local character. As a consequence, the 
sequential updating of the Trotter time slices, which we 
mentioned in Sec. El can no longer be exploited to reduce 
the numerical effort for the evaluation of the fermionic 
weight. However, in combination with the reweighting 
method, the latter is only calculated when measurements 
of observables are made. The enormous advantage of the 
PCs, leading to completely uncorrelated phonon config- 
urations, clearly outweighs this drawback. Nevertheless, 
this restriction has to be kept in mind when consider- 
ing possible extensions to many-electron systems. Apart 
from this, the PC representation can also be applied to 
the more general case of more than one electron, since 
the bosonic action [Eq. (|46[) . on which the transforma- 
tion relies, remains unchanged]. This even holds for the 
case of more general models including, e.g., spin-spin or 
Hubbard-type interactions, as long as the phonon oper- 
ators enter in the same form as in the Holstein model. 

Another important point is the combination of the 
PCs with the reweighting method. Using the latter, the 
changes to the original momenta p, which are made in the 
simulation, do not depend in any way on the electronic 
degrees of freedom. Thus we are actually sampling a set 
of N independent harmonic oscillators, as described by 
the purely bosonic action Sh- The crucial requirement 
for the success of this method is the use of the Lang- 
Firsov transformed model, in which the polaron effects 
are separated from the zero-point motion of the oscilla- 
tors around their current equilibrium positions. 

Finally, for the untransformed model, Eq. Q), the 
bosonic action can be obtained from Eq. I|33[) by replac- 
ing p with x (see Sec. IV All and a transformation to PCs 
could also be used. However, as discussed in Sec. I VII 
without the Lang-Firsov transformation, the reweighting 
procedure fails. Consequently, using the standard ap- 
proach, the phonon coordinates x would depend on the 
electronic degrees of freedom, and this makes exact sam- 
pling impossible for the untransformed model. 



VIII. RESULTS: VPA 

In order to test the validity of VPA we calculated 
the total energy [Eq. i|18|) ] and the quasiparticle weight 
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N = 4 




FIG. 2: Total energy E (top) and quasiparticle weight zo 
(bottom) as functions of the electron-phonon coupling A for 
different values of the phonon frequency u>. Symbols corre- 
spond to VPA results, while full lines represent exact T — 
data obtained with the Lanczos method fRef. FFfTl . Dashed 
lines are results of the HLF approximation. 




FIG. 3: Total energy E as a function of the electron-phonon 
coupling g (see text) for different values of the phonon fre- 
quency u>. Symbols correspond to VPA results, while full 
lines represent data obtained with the globallLocal method 
(Ref. l57l) . The dashed line represents the atomic- limit result 
(p — oo). 
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FIG. 4: Polaron-size parameter 7,5 as a function of the 
electron-phonon coupling A for various distances 8. The pa- 
rameter 7 of the standard Lang-Firsov transformation (see 
Sec. II11I is also shown. 

[Eq. on a cluster of four sites for various phonon fre- 
quencies w and compared the results with those of Mar- 
siglio obtained by Lanczos diagonalizationiSi The com- 
parison is depicted in Fig. [21 The values of ui have been 
chosen to lie in the nonadiabatic regime to > 1 where the 
zero-phonon approximation of the VPA is sensible. The 
overall agreement is strikingly good. Minor deviations 
from the exact results increase with decreasing phonon 
frequency. For the smallest frequency shown, Q = 1.0, 
the curve for the HLF approximation is also depicted. 
It reveals that VPA represents a significant improvement 
over the HLF approximation, underlining the importance 
of the extended polaron cloud. 

The comparison with exact results obtained with Lanc- 
zos was restricted to small clusters with N = 4 in order 
to achieve convergence with respect to the number of 
phonon states included in the calculation (see Sec.[H]). To 
further scrutinize the accuracy of the VPA we also com- 
pare the results of the latter for the total energy with the 
variational global-local method which has been shown to 
give reliable results over a large range of parameters^ 
We chose N — 32 for which finite-size effects are already 
very small (see Sec. IIX|I. Moreover, following Romero et 



13 



al.^L in Fig. we plot E/u) over g with g = y / XW/(2uj). 
Similar to the case N = 4 shown in Fig.|21we find a very 
good agreement for large values of uj over the whole range 
of electron-phonon coupling, whereas for smaller ui the 
VPA results begin to bend away from the correct curve 
and collapse to the strong-coupling, atomic-limit result 
for large g. We would like to point out that the maximum 
electron-phonon coupling strength in Fig. [5] corresponds 
to A w 40 (for uj = 4.0), in contrast to Fig. [5] where 
A < 2. Figures El and reveal that in the nonadiabatic 
regime uj 1 VPA yields a very good agreement with 
the exact data and the Global-Local method even in the 
intermediate and strong-coupling regime. This behavior 
can easily be understood considering the assumptions of 
the VPA. The zero-phonon approximation becomes ex- 
act in the nonadiabatic limit Q — > oo, where the energies 
of phonon excitations are too high to have an effect on 
the ground state. Finally, we would like to mention the 
possibility of comparing the VPA with the QMC results 
presented in the following section. This has been done 
for a variety of parameters, but we have found that it 
is difficult to distinguish between deviations due to the 
shortcomings of the VPA and due to temperature effects 
in the QMC results. Consequently, we have decided to 
confront the VPA with another approved ground-state 
method, namely, the global-local method, which gives a 
much clearer picture. 

In Fig. 2] we show results for the variational displace- 
ment fields jg, which give us a measure for the size of 
the polaron. For uj = 0.1 we see an abrupt crossover 
from a large to a small polaron at A w 1.2. For smaller 
values of the coupling, the electron induces lattice dis- 
tortions at neighboring sites even at a distance of more 
than three lattice constants. Above Aw 1.2 we have a 
mobile small polaron extending over a single site only. 
In contrast, for a larger value of the phonon frequency 
uj = 4.0, there is no crossover and we have a somewhat 
extended (large) polaron even for large values of A. The 
same behavior has been found by Marsiglio 3 ^ who deter- 
mined the correlation function (mxi+s) by Lanczos diag- 
onalization for a restricted phonon basis. Within VPA 
we have the relation (niXj+s) = 75. The main difference 
is that in Marsiglio's results, the crossover to a small 
polaron for Q = 0.1 occurs at a smaller value of the cou- 
pling Awl. Nevertheless, the simple VPA reproduces 
the main characteristics of the transition of the Holstcin 
polaron as the coupling strength A is increased. Finally 
Fig. Q] also shows the result for the parameter 7 of the 
standard Lang-Firsov transformation (see Sec. IIII|) . For 
Q = 0.1, the curves for 7 and 75=0 are identical above 
the critical value A ~ 1.2. This is not surprising since, 
in this regime, we have a small polaron extending over 
a single site only, which is well described by the local 
Lang-Firsov transformation defined in Eq. (J5J . For larger 
values of the phonon frequency (see Fig.^J, 7 and 70 do 
not coincide above a critical value of the coupling, but 
the difference vanishes as A — > 00 . In contrast to the adi- 
abatic regime, the polaron remains an extended object 



up to very large values of the coupling, so that the local 
ansatz of the Lang-Firsov transformation does then not 
provide the correct description for finite values of A (see 
also Ref.il. 

IX. RESULTS: QUANTUM MONTE CARLO 

As our approach is based on a discretized imaginary 
time, it is important to study the convergence of any 
results with increasing number of time slices, L, which 
determines the error due to the Suzuki- Trotter approxi- 
mation of Eq. (|25(l . L was chosen such that systematic 
errors are smaller than the statistical errors of the re- 
sults. For all observables considered here we have found 
the usual (At) 2 dependence of the Suzuki- Trotter error. 
Depending on the phonon frequency u> we have found 
values of At = 1/30 (for u < 1) and At = 1/40 (for 
uj > 1) to be sufficient even for the most accurate results 
of this paper. Moreover, as indicated in Fig. EI error bars 
for the QMC data presented are always smaller than the 
symbols used in the figures and are therefore not shown. 
Finally, lines connecting data points obtained with QMC 
in Figs. I^ITTTlare guides to the eye only. 

To test our QMC algorithm we have performed several 
comparisons with other methods. First, we have checked 
that the QMC reproduces the exact results obtained with 
Lanczos on a four-site cluster. Apart from temperature 
effects, an excellent agreement has been found for several 
different values of the phonon frequency. Second, as the 
QMC results are all for finite temperature, we have also 
compared them with an exact solution for the two-site 
system, which is valid for arbitrary temperature. We 
have found a perfect agreement over the whole range of 
values for /?, uj, and A, and can therefore exclude the 
possibility of any systematic errors. 

A. Kinetic energy 

We begin our discussion of the results with the kinetic 
energy of the electron, given by Eq. (|38|l . which has previ- 
ously been calculated by several authorsi 22 i 23 i 25 i 43 i 49 i 57 i 65 
In Fig. we show results for E^ on a 32-site cluster, with 
[3t = 10 and for several values of the phonon frequency. 
While for small values of Q there is a rapid decrease of 
the absolute value of the kinetic energy in the vicinity of 
A = 1, the cross over becomes smoother as uj increases. 
This agrees with the findings of previous studies and re- 
sembles closely the behavior of the total energy discussed 
above. For large values of A and Q < 1 we find E^ ~ A" 1 
as predicted by small-polaron theory. 71 This contrasts 
strongly with the behavior of the quasiparticle weight zq 
[see Fig. EJb)] which decreases much faster and is ex- 
ponentially suppressed in the small-polaron regime^ As 
pointed out by Fehske et al.^ for the case of the Holstein 
model, the quasiparticle weight is exactly the inverse of 
the ratio m e ff/m where m e ff and m denote the effective 
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FIG. 5: Negative kinetic energy as a function of the 
electron-phonon coupling A for various values of the phonon 
frequency Q. 
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FIG. 6: Negative kinetic energy as a function of the electron- 
phonon coupling A for various values of the inverse tempera- 
ture j3 for various values of the phonon frequency Q. 



and free mass of the electron, respectively. Hence, in the 
small-polaron regime, the effective mass increases expo- 
nentially, while the kinetic energy still has a finite value. 
We ascribe this behavior to the undirected motion of the 
electron inside the phonon cloud, which gives rise to a 
nonzero kinetic energy even for large values of A. How- 
ever, since the polaron bandwidth is exponentially nar- 
rowed with increasing A, the polaron is almost localized. 

To study the influence of temperature we have calcu- 
lated the kinetic energy for a system of 32 sites, with 
Q = 1.0 and for various values of fit (see Fig. EJJ. As fit 
increases, \E^\ increases for A < 2. However temperature 
effects are obviously very small in the strong-coupling 
regime. For fit = 1, \Ek\ decays in a smooth way as A is 
increased, while for lower temperatures we find the typi- 
cal rather abrupt crossover near A = 1, as in Fig. [5j De 
Raedt and LagendijkH have calculated the kinetic energy 
for the same set of parameters using their QMC method. 
However, the lowest temperature they could reach was 
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FIG. 7: Negative kinetic energy as a function of the electron- 
phonon coupling A for different numbers of lattice sites N. 



1.6 



1.4 



1.2 



, . 1 1 1 
(a) 


iii 

(O=l.0,X=l 




o-^N = 6 




^N= 8 




b-hN=16 




c^3N = 32 


i i 


iii 



0.1 0.2 0.3 0.4 

2r 



0.5 0.6 
l/P 



0.7 0.8 0.9 1 



1.5 



pq 



'(b) 



co= 1.0, X= 1 

0-o 0t= 1 

B-ElPt= 3 

o~o Pt= 5 
^Pt= 10 




0.03125 



0.0625 



9 9 



0.1 0.125 



0.25 



1/N 
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the number of lattice sites N. 



fit = 5 which, according to Fig. H3 is still quite different 
from the ground-state result. Moreover, their calcula- 
tions did not include dynamical effects of the phonon de- 
grees of freedom. As a consequence, for fit = 1, they do 
not obtain the correct behavior of the kinetic energy as 
a function of A. Finally, Romero et al& and Jeckelmann 
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and White 49 calculated the kinetic energy for T = on 
a 32-site cluster and for an infinite system, respectively. 
Their results are in a good agreement with our findings, 
although small deviations due to temperature and finite- 
size effects are visible. Nevertheless, we can conclude 
from Fig. El that a value of (3t = 10 should be sufficient 
to obtain results which are representative of the ground 
state. 

We now turn our attention to finite-size effects. In 
Fig. 0we show the kinetic energy for uj = 0.5, (3t = 10, 
and for various number of lattice sites. For N > 16 the 
results for E^ are well converged over the whole range of A 
and finite-size effects are very small. This agrees with the 
findings of other authors i^ 3 * 3 ^^ Figure |H1 shows the ki- 
netic energy as a function of temperature, for Cu — 1.0 and 
various numbers of lattice sites N. Moreover we chose 
A = 1 , as the influence of the system size is largest in the 
cross over regime. Figure |HJa) clearly demonstrates that 
finite-size effects are most pronounced at low tempera- 
tures, while they are completely smeared out at higher 
temperatures, since high-temperature properties are de- 
termined by integral quantities, such as energy moments 
(E v ), which have a small size dependence, while low- 
temperature features are governed by energetically low- 
lying eigenvectors. To further illustrate the influence of 
the system size, we plot in Fig. Efb) the negative kinetic 
energy as a function of 1/N again for w = A = 1 and 
for various values of (3. As before, error bars are smaller 
than the symbol size, but due to the very high accu- 
racy of the data, the systematic errors due to the finite 
number of Trotter slices L are comparable to the statisti- 
cal errors. The results show that very good convergence 
with respect to the number of lattice sites is achieved 
for rather small N. In fact, for the highest tempera- 
ture shown ((3t =1), the line connecting the data points 
becomes vertical already at N — 8, while for (3t = 10 
convergence is reached for N = 16. Hence, if we con- 
sider these findings in the context of the usual finite-size 
scaling analysis where one plots the data as a function of 
a suitably chosen power of 1/N and extrapolates to the 
infinite system (i.e., 1/N — > 0), we have here the spe- 
cial case of a linear dependence with zero slope at large 
enough N. Thus, in contrast to the half-filled Holstein 
model of spinless fermion, for which a finite-size anal- 
ysis has been performed by two groups^*^ we merely 
find that the results converge within the accuracy of our 
calculations already for rather small systems. 



B. Total energy 

Next we consider the total energy E, given by Eq. i|4ip. 
In Fig. |5| we present the total energy for a cluster of 
32 sites and various values of the phonon frequency. 
Finite temperature effects increase as we approach the 
low- frequency regime, and for uj = 0.1 we clearly see a 
strong deviation from the ground-state result E = — 2t 
for A = 0. The frequency-dependence of the temperature 




FIG. 9: QMC results for the total energy E as a function 
of the electron-phonon coupling A for various values of the 
phonon frequency Q. Here and in subsequent figures lines are 
guides to the eye, and errorbars for the QMC data are smaller 
than the symbols shown. 

effects can easily be understood if we consider the exact 
result for the kinetic energy of N independent harmonic 
oscillators 

uj 2 . Nu (1 1 \ 
^ph=2£<P?> = — {-2 + ^—,) ■ ( 47 ) 

which is identical to the second term in Eq. (|41[) . For low 
temperatures we have (p 2 ) » 0.5 + e _/3w , with a correc- 
tion to the ground-state value of 0.5 that increases with 
decreasing ui. These temperature effects on E due to the 
oscillator energy do not depend on A [see Eq. (|T7|) ] and 
therefore shift the total energy curve by the same amount 
for all values of the coupling. A comparison with the dis- 
cussion of the kinetic energy reveals that temperature 
effects are much smaller for other observables due to the 
absence of the strongly temperature-dependent phonon 
energy terms P p and E' ph [see Eqs. Q and l|4"TJl ]. 

The dependence on u> agrees well with existing 
work, 22 i 23 . 25 . 31 i 36 i 38 i 41 i 47 i 56 i 57 i 59 i 70 It is known^ that at 
zero temperature and for small values of the phonon 
frequency, u> < 0.5, the total energy displays a rather 
sharp transition around A » 1, where the cross over 
from a large to a small polaron occurs. In ED studies 
of small clusters— a kink in E has been observed, which 
is smeared out in the finite-temperature QMC results. 
Nevertheless, we observe the same rounding of the en- 
ergy curve with increasing uj (Ref. I36T) . As discussed by 
Marsiglic 36 the kink in the total energy is merely a finite- 
size effect. As the system size increases the discontinuity 
disappears, in accordance with the fact that the ground 
state of the Holstein polaron is an analytic function of 
the coupling parameter A (Ref. 0). 

Finally, it is interesting to note that in contrast to the 
kinetic energy — E^, which shows a sharp decrease near 
A = 1 in the adiabatic regime (see, e.g., Fig. [HJ, the to- 
tal energy does not change significantly. As discussed for 
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the two-dimensional case by Kornilovitch^ this can be 
explained as follows. For small ratios to/t, the phonon 
energy associated with the term P of Hamiltonian (JJJ 
is small and the system is governed by the balance of 
the electronic kinetic energy and the energy due to the 
electron-phonon coupling. In the transformed model, the 
latter is given by E p as defined by Eq. JSJ . When the ra- 
tio of the two energies approaches unity (equivalent to 
A = 1), it becomes energetically favorable for the elec- 
tron to localize (losing kinetic energy) and increase its 
potential energy. This leads to finite displacements of 
the oscillators in the vicinity of the electron and increases 
the potential energy of the phonons. Hence, near A = 1 
the energy of the system is redistributed from kinetic to 
potential energy so that E remains almost unchanged. 
This is exactly what we see in Fig. 
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FIG. 10: Momentum distribution n(k) as a function of A for 
various wave vectors k. 



C. Momentum distribution and oscillator momenta 



Following Zhang et alw* we also calculated the mo- 
mentum distribution n(£:), given by Eq. (|42|) . for differ- 
ent wave vectors k (Fig. llOjl . To compare with their 
DMRG 7 ? results we chose the same parameters N = 6 
and ui = 1.0. Moreover, we took j3t = 10 since the cal- 
culations of Zhang et al. were for the ground state. For 
A = the ground state has momentum k = 0, so we 
have n(0) = 1 and n(k ^ 0) = 0. With increasing 
coupling n(0) decreases in a way similar to the kinetic 
energy (cf. Fig. [5}, while n(k) for k ^ increases. In 
the strong-coupling limit A — > oo, n(k) approaches the 
value l/N =1/6 for all k. This is a simple consequence 
of the localization of the electron for A = oo. Although 
the curve for k — looks very similar to the results of 
Zhang et al. we find a slightly stronger decrease of n(0) 
in the intermediate coupling regime. This deviation is 
no temperature effect of our QMC method but proba- 
bly originates from the fact that Zhang et al. obtained 
their results for n(0) by integrating over an approximate 
spectral function. 

In Sec. IVll we mentioned that, within the Lang-Firsov 
approach, the phonon degrees of freedom only show a 
weak dependence on the electron-phonon coupling, in 
contrast to the standard approach, where the average 
oscillator coordinate (x) increases strongly with A due to 
the displacement in the presence of an electron. The 
weak dependence of the vibrational energy of the lo- 
cal oscillator, which is proportional to (p 2 ), on A is 
shown in Fig. ^] For A = we have the result (p 2 ) = 
0.5 + [exp(/3u>) — 1] _1 [see Eq. (|47f) ] for a free oscillator. In 
the intermediate coupling regime, (p 2 ) takes on a mini- 
mum, corresponding to a reduction of merely 4% and ap- 
proaches the value for A = again in the strong-coupling 
limit. As the Lang-Firsov transformation does not af- 
fect the phonon momenta p (see Sec. Illlf) . the result for 
(p 2 ) as a function of A is the same in the untransformed 
Holstein model. However, the significant advantage of 
the proposed method is that the phonon momenta are 




FIG. 11: Mean square of the phonon momentum p as a 
function of the electron-phonon coupling A. 



sampled instead of the coordinates x. Thus the proba- 
bility distribution associated with the degrees of freedom 
to be sampled has only a small variance compared to 
the standard method, which makes the simulations much 
more effective. The dependence of (p 2 ) on the coupling 
strength A and the temperature has first been studied by 
Ranninger and Thibblinii for the two-site polaron prob- 
lem. For such a small system, the minimum of (p 2 ) is even 
more pronounced, while for larger systems the average ef- 
fect of the electron on a local oscillator is more and more 
washed out. Ranninger and Thibblin^L ascribed the de- 
viation of the vibrational energy from the free-oscillator 
result to anharmonic effects, which are visible only at low 
enough temperatures. This can clearly be seen in Fig. 1111 
where the minimum of (p 2 ) becomes less pronounced and 
is shifted to smaller values of A as the temperature in- 
creases. 
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D. Performance 

We conclude this section with a discussion of the per- 
formance of the QMC approach. From the results pre- 
sented above it is obvious that the method enables us to 
study a very wide range of parameters. Hence, for exam- 
ple, we have performed calculations for 0.1 < Q < 4.0 (see 
Fig. 0) . Simulations in the adiabatic regime would be ex- 
tremely difficult within the standard approach, since the 
autocorrelation times grow as (wAt) -2 . However, in ma- 
terials such as the manganites, the frequencies of the rel- 
evant phonon modes are known to be small (<D < 0.5, see, 
e.g., Ref.0) so that our method could represent an im- 
portant step forward towards the simulation of electron- 
phonon models with realistic parameters. Also, we are 
able to reach very low temperatures (it < 20 and clus- 
ters large enough to avoid finite-size effects with modest 
computational effort. Another key advantage is that the 
method becomes more and more efficient as the coupling 
strength A increases, which is due to the use of the Lang- 
Firsov transformation. In our results we find that statis- 
tical errors of expectation values of phonon operators are 
larger than, e.g., the errors of the kinetic energy. Finally, 
the errors increase slightly as we approach the adiabatic 
and/or low-temperature regime Co — > and fit — > oo, 
respectively. 

To demonstrate the efficiency of our method we give 
some figures for the CPU time of the simulations. A 
typical QMC run for 32 lattice sites, fit = 5, Q = 1.0, 
and A ~ 1 (i.e., near the small-polaron crossover) only 
takes 5 min of CPU time on a 650 MHz Pentium III 
PC. For such a run relative errors of, for example, the ki- 
netic energy are less than 1.0%. Away from the crossover 
point, the same accuracy can be obtained within a few 
seconds. For fit = 10, the temperature used in most 
of the calculations presented in this paper, a QMC run 
with A near the crossover value and with similar statisti- 
cal errors as mentioned above takes about 80 min of CPU 
time. Hence, although not as efficient as the specialized 
one-electron methods; 2 ^ 2 ^* 2 ^^* 3 ^ 3 ^ our approach signifi- 
cantly reduces the numerical effort compared to previous 
methods which were often run on supercomputers and 
did not reach the parameters (low temperature and small 
phonon frequency) and accuracy of our simulations. 

X. CONCLUSIONS 

We have presented a simple variational approach to the 
Holstein model, which incorporates an extended Lang- 
Firsov transformation. This approach is easily applicable 
to infinite systems and represents a marked improvement 
over the standard small-polaron approximation, which is 
only useful in the nonadiabatic, strong-coupling regime. 

More importantly, we have introduced an exact QMC 
method for the Holstein model, which is based on the 
standard Lang-Firsov transformation of the Hamiltonian. 
The phonon momenta are represented in terms of princi- 



pal components, which enables us to sample completely 
uncorrelated configurations, while the electronic degrees 
of freedom are taken into account exactly by use of a 
reweighting method for calculating observables. Thereby, 
we avoid the numerically expensive evaluation of the elec- 
tronic weights in the updating process. The present ap- 
proach can be applied for a wide range of parameters 
with relatively small computational effort. In particular, 
efficient simulations can be performed in the adiabatic 
regime cD < 1, which is of special interest in connec- 
tion with materials such as the manganites. In the one- 
dimensional case considered here, a sign problem result- 
ing from the Lang-Firsov transformation on small sys- 
tems has been found to have only a small effect on the 
statistics. Tests have been presented in the one-electron 
case and reveal that the method reproduces Lanczos di- 
agonalization results in the regime where the latter are 
applicable, namely, for very small systems, small to mod- 
erate electron-phonon coupling and for sufficiently large 
phonon frequency. Moreover, a satisfactory agreement 
with other methods has been found. Owing to the exact 
treatment of the electronic degrees of freedom and the 
sampling of the phonons, the method is free of any auto- 
correlations. The use of the Lang-Firsov transformation, 
which is essential for the applicability of the reweighting 
method, substantially improves the statistics, allowing 
for very accurate results. 

Despite the great computational efficiency of our 
method compared to the standard approach, even faster 
methods exist. For example, the QMC simulations of dc 
Raedt and Lagendijk 22 i 23 i 24 and Kornilovitch 25 ! 31 1 32 seem 
to be numerically faster due to the analytic integration 
over the phonon degrees of freedom which significantly 
reduces statistical errors. However, both methods are re- 
stricted in their applicability as discussed in Sec. [I] In 
particular, an extension to many-electron systems seems 
impossible, since simulations will be restricted by a severe 
minus sign problem similar to other world-line methods. 
In contrast, the method presented here is not restricted 
to the single-electron limit in principle, although some 
modifications will be necessary. As pointed out in previ- 
ous sections, most of the ideas proposed here, such as the 
use of the transformed model, the reweighting method, 
and the PC representation, remain unchanged if we con- 
sider more than one electron. The required modifications 
concern mainly the fermionic weight Wf [Eq. I|29l) ]. There 
is a Hubbard-like interaction term coming from the Lang- 
Firsov transformation (see Sec. 1111(1 . and the one-electron 
basis states used here (Eq. IKSH have to be replaced by 
the corresponding set of many-electron states. Since the 
number of such basis states, and thereby the linear di- 
mension of the matrices fi, n, and D (see Sec. IV A|) . grows 
exponentially with the system size, an exact treatment 
of the fermion degrees of freedom will become increas- 
ingly difficult. Consequently, a more refined approach 
based on, e.g., the use of determinant methods will be re- 
quired. For the bipolaron problem of two electrons with 
opposite spin, on which work is currently in progress, 
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the computational effort can be significantly reduced by 
exploiting the conservation of the total quasimomentum 
(see, e.g., Ref.Esh. leading to a computer time that grows 
with the cube of the system size. The more general case 
of, e.g. , a quarter- filled band corresponding to the colos- 
sal magnetoresistance regime of the manganites, requires 
further consideration, and the effect of the sign problem 
remains to be investigated. Moreover, the performance of 
such an approach has to be compared to existing many- 
electron QMC methods for the Holstein model. Finally, 
the method can be generalized to more complicated mod- 
els including, e.g., a coupling of the electrons to local 
spins as in the Kondo or double-exchange model for the 
manganites. 



Acknowledgments 



This work was partially supported by the Austrian Sci- 
ence Fund (FWF), project No. P15834. M.H. is grateful 
to the Austrian Academy of Sciences for financial sup- 
port. We would like to acknowledge helpful discussions 
with Markus Aichhorn, Holger Fehske, Winfried Koller, 
and Alexander Priill. We would also like to thank Frank 
Marsiglio and Aldo Romero for providing us with some 
of the data presented in this paper. 



* Electronic address: hohenadler@itp.tu-graz.ac. at 

1 Y. Bar- Yam, J. Mustre de Leon, and A. R. Bishop, 
eds., Lattice Effects in High Temperature Superconductors 
(World Scientific, Singapore, 1992). 

2 D. M. Edwards, Adv. Phys. 51, 1259 (2002). 

3 A. J. Millis, R. Mueller, and B. I. Shraiman, Phys. Rev. B 
54, 5405 (1996). 

4 E. Dagotto, S. Yunoki, and A. Moreo, in Physics of Man- 
ganites, edited by T. A. Kaplan and S. D. Mahanti (Kluwcr 
Publishing, 1999), p. 39. 

5 T. Holstein, Ann. Phys. (N.Y.) 8, 325; 8, 343 (1959). 

6 H. Froehlich, H. Pelzer, and S. Zienau, Philos. Mag. 41, 
221 (1950). 

7 W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. 
Lett. 42, 1698 (1979). 

8 R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys- 
ical Review D 24, 2278 (1981). 

9 D. J. Scalapino and R. L. Sugar, Phys. Rev. B 24, 4295 

(1981) . 

10 G. Levine and W. P. Su, Phys. Rev. B 42, 4143 (1990). 

11 G. Levine and W. P. Su, Phys. Rev. B 43, 10 413 (1991). 

12 P. Niyaz, J. E. Gubernatis, R. T. Scalettar, and C. Y. 
Fong, Phys. Rev. B 48, 16 011 (1993). 

13 J. E. Hirsch, D. J. Scalapino, R. L. Sugar, and R. Blanken- 
becler, Phys. Rev. Lett. 47, 1628 (1981). 

14 J. E. Hirsch, R. L. Sugar, D. J. Scalapino, and R. Blanken- 
becler, Phys. Rev. B 26, 5033 (1982). 

15 J. E. Hirsch and E. Fradkin, Phys. Rev. Lett. 49, 402 

(1982) . 

16 E. Fradkin and J. E. Hirsch, Phys. Rev. B 27, 1680 (1983). 

17 J. E. Hirsch and E. Fradkin, Phys. Rev. B 27, 4302 (1983). 

18 W. von der Linden, Phys. Rep. 220, 53 (1992). 

19 R. T. Scalettar, N. E. Bickers, and D. J. Scalapino, Phys. 
Rev. B 40, 197 (1989). 

20 R. M. Noack, D. J. Scalapino, and R. T. Scalettar, Phys. 
Rev. Lett. 66, 778 (1991). 

21 F. Marsiglio, Phys. Rev. B 42, 2416 (1990). 

22 H. De Raedt and A. Lagendijk, Phys. Rev. Lett. 49, 1522 

(1982) . 

23 H. De Raedt and A. Lagendijk, Phys. Rev. B 27, 6097 

(1983) . 

24 H. De Raedt and A. Lagendijk, Phys. Rev. B 30, 1671 

(1984) . 

25 P. Kornilovitch, J. Phys.: Condens. Matter 9, 10675 



(1997). 

26 R. P. Feynman, Physical Review 97, 660 (1955). 

27 H. de Raedt and A. Lagendijk, Z. Phys. B: Condens. Mat- 
ter 65, 43 (1986). 

28 P. E. Kornilovitch and E. R. Pike, Phys. Rev. B 55, R8634 

(1997) . 

29 N. V. Prokof'ev and B. V. Svistunov, Phys. Rev. Lett. 81, 
2514 (1998). 

30 A. S. Mishchenko, N. V. Prokof'ev, A. Sakamoto, and B. V. 
Svistunov, Phys. Rev. B 62, 6317 (2000). 

31 P. E. Kornilovitch, Phys. Rev. Lett. 81, 5382 (1998). 

32 P. E. Kornilovitch, Phys. Rev. B 60, 3237 (1999). 

33 E. Berger, P. Valasek, and W. von der Linden, Phys. Rev. 
B 52, 4806 (1995). 

34 R. H. McKenzie, C. J. Hamer, and D. W. Murray, Phys. 
Rev. B 53, 9676 (1996). 

35 P. Sengupta, A. W. Sandvik, and D. K. Campbell, Phys. 
Rev. B 67, 245103 (2003). 

36 F. Marsiglio, Physica C 244, 21 (1995). 

37 V. V. Kabanov and D. K. Ray, Phys. Lett. A 186, 438 
(1994). 

38 A. S. Alexandrov, V. V. Kabanov, and D. K. Ray, Phys. 
Rev. B 49, 9915 (1994). 

39 I. G. Lang and Y. A. Firsov, Sov. Phys. JETP 16, 1301 
(1962). 

40 G. D. Mahan, Many-particle Physics (Plenum Press, New 
York, 1990), 2nd ed. 

41 J. Ranninger and U. Thibblin, Phys. Rev. B 45, 7730 
(1992). 

42 F. Marsiglio, Physics Letters A 180, 280 (1993). 

43 G. Wellein, H. Roder, and H. Fehske, Phys. Rev. B 53, 
9666 (1996). 

44 G. Wellein and H. Fehske, Phys. Rev. B 56, 4513 (1997). 

45 J. M. Robin, Phys. Rev. B 56, 13 634 (1997). 

46 M. Capone, W. Stephan, and M. Grilli, Phys. Rev. B 56, 
4484 (1997). 

47 M. Capone, S. Ciuchi, and C. Grimaldi, Europhys. Lett. 
42, 523 (1998). 

48 E. Jeckelmann and S. R. White, Phys. Rev. B 57, 6376 

(1998) . 

49 E. Jeckelmann, C. Zhang, and S. R. White, Phys. Rev. B 
60, 7950 (1999). 

50 C. Zhang, E. Jeckelmann, and S. R. White, Phys. Rev. B 
60, 14 092 (1999). 



19 



51 A. Weifie, H. Fehske, G. Wellein, and A. R. Bishop, Phys. 
Rev. B 62, R747 (2000). 



52 
53 



W. Stephan, Phys. Rev. B 54, 8981 (1996). 
M. Hohenadler, M. Aichhorn, and W. von der Linden, 
cond-mat/0306740 (2003), unpublished. 

54 J. Bonca, S. A. Trugman, and I. Batistic, Phys. Rev. B 60, 
1633 (1999). 

55 L. C. Ku, S. A. Trugman, and J. Bonca, Phys. Rev. B 65, 
174306 (2002). 

56 A. H. Romero, D. W. Brown, and K. Lindenberg, Phys. 
Rev. B 59, 13 728 (1999). 

57 A. H. Romero, D. W. Brown, and K. Lindenberg, Phys. 
Rev. B 60, 4618 (1999). 

58 V. Cataudella, G. De Filippis, and G. Iadonisi, Phys. Rev. 
B 63, 52 406 (2001). 

59 H. Fehske, J. Loos, and G. Wellein, Phys. Rev. B 61, 8016 
(2000). 

60 O. S. Barisic, Phys. Rev. B 65, 144301 (2002). 

61 O. S. Barisic, cond-mat/0211189 (2002), unpublished. 

62 M. Acquarone, M. Cuoco, C. Noce, and A. Romano, Phys. 
Rev. B 63, 035110 (2001). 

63 S. Ciuchi, F. de Pasquale, S. Fratini, and D. Feinberg, 
Phys. Rev. B 56, 4494 (1997). 



64 H. Lowen, Phys. Rev. B 37, 8661 (1988). 

65 E. V. L. de Mello and J. Ranninger, Phys. Rev. B 55, 14 
872 (1997). 

66 G. G. Batrouni and R. T. Scalettar, in Quantum Monte 
Carlo Methods in Physics and Chemistry, edited by M. P. 
Nightingale and C. J. Umrigar (Kluwer Academic Publish- 
ers, 1998), p. 65. 

67 J. N. Kapur and H. K. Kesavan, Entropy Optimization 
Principles with Applications (Academic Press, 1992). 

68 A. C. Davison and D. V. Hinkley, Bootstrap Methods and 
their Application (Cambridge University Press, 1997). 

69 W. H. Press, Numerical recipes in Fortran 77, 
http://www.numrec.com. 

7(1 F. Marsiglio, in Recent Progress in Many-Body Theories, 
edited by E. Schachinger, H. Mitter, and H. Sormann 
(Plenum Press, 1995), vol. 4. 

71 A. S. Alexandrov and N. F. Mott, Rep. Prog. Phys. 57, 
1197 (1994). 

72 A. Weifie and H. Fehske, Phys. Rev. B 58, 13 526 (1998). 

73 The coupling parameter g used in Ref.[H(3is related to ours 
hy \ = 2g 2 /{ujt). 



